function [Nbot] = TEC_BOT(h, A1, A2, A3, hmF2, hmF1, hmE, B2bot, B1top,...
                          B1bot, BEtop, BEbot)
%
% Bottomside Electron Density
%
%DESCRIPTION:
%This function computes the Bottomside Electron Density.
%
%PROTOTYPE:
% [Nbot] = TEC_BOT(h, A1, A2, A3, hmF2, hmF1, hmE, B2bot, B1top, B1bot, BEtop, BEbot)
%
%--------------------------------------------------------------------------
% INPUTS:
%   h          [1x1]       Height                    [km]
%   A1         [1x1]       F2-Layer Amplitude        [1e11 m-3]
%   A2         [1x1]       F1-Layer Amplitude        [1e11 m-3]
%   A3         [1x1]       E-Layer Amplitude         [1e11 m-3]
%   hmF2       [1x1]       F2-Layer Max. Den. Height [km]
%   NmF1       [1x1]       F1-Layer Max. Density     [1e11 m-3]
%   NmE        [1x1]       E-Layer Max. Density      [1e11 m-3]
%   hmF1       [1x1]       F1-Layer Max. Den. Height [km]
%   hmE        [1x1]       E-Layer Max. Den. Height  [km]
%   B2bot      [1x1]       F2-Layer Bottom Thickness [km]
%   B1top      [1x1]       F1-Layer Top Thickness    [km]
%   B1bot      [1x1]       F1-Layer Bottom Thickness [km]
%   BEtop      [1x1]       E-Layer Top Thickness     [km]
%   BEbot      [1x1]       E-Layer Bottom Thickness  [km]
%--------------------------------------------------------------------------
% OUTPUTS:
%   Nbot       [1x1]       Bottomside Electron Dens. [m-3]
%--------------------------------------------------------------------------
%
%NOTES:
% (none)
%
%CALLED FUNCTIONS:
% (none)
%
%UPDATES:
% (none)
%
%REFERENCES:
% [1] "Ionospheric Correction Algorithm for Galileo Single-Frequency Users"
%      - European GNSS (Galileo) Open Service
% [2] "Electron Density Models and Data for Transionospheric Radio
%      Propagation" - Report ITU-R P.2297-1 (05/2019)
%
%AUTHOR(s):
%Luigi De Maria, Matteo D'Addazio, 2022
%

%% Main Code

%Pre-Settings
A = [A1; A2; A3];

%Height Case Determination
if     h > hmE
    BE = BEtop;
elseif h <= hmE
    BE = BEbot;
end
if     h > hmF1
    BF1 = B1top;
elseif h <= hmF1
    BF1 = B1bot;
end

%Exponential Argument for Each Layer
alpha = zeros(3,1);
%Initial Check
if h >= 100 %[2] says >= 90 km
    alpha(1) = (h - hmF2) / (B2bot);
    alpha(2) = (h - hmF1) / (BF1) * exp(10/(1 + abs(h - hmF2)));
    alpha(3) = (h - hmE)  / (BE)  * exp(10/(1 + abs(h - hmF2)));
else
    alpha(1) = (100 - hmF2) / (B2bot);
    alpha(2) = (100 - hmF1) / (BF1) * exp(10/(1 + abs(100 - hmF2)));
    alpha(3) = (100 - hmE)  / (BE)  * exp(10/(1 + abs(100 - hmF2)));
end

%Electron Density
s = zeros(3,1);
for i = 1 : 3
    if     abs(alpha(i)) > 25
        s(i) = 0;
    elseif abs(alpha(i)) <= 25
        s(i) = A(i) * (exp(alpha(i)))/(1 + exp(alpha(i)))^2;
    end
end

if     h >= 100
    Nbot = (s(1) + s(2) + s(3)) * 1e11;
elseif h < 100
    %Cerrective Terms
    %Term 1
    if     abs(alpha(1)) > 25
        ds1 = 0;
    elseif abs(alpha(1)) <= 25
        ds1 = 1/B2bot * (1 - exp(alpha(1))/(1 + exp(alpha(1))));
    end
    %Term 2
    if     abs(alpha(2)) > 25
        ds2 = 0;
    elseif abs(alpha(2)) <= 25
        ds2 = 1/BF1 * (1 - exp(alpha(2))/(1 + exp(alpha(2))));
    end
    %Term 3
    if     abs(alpha(3)) > 25
        ds3 = 0;
    elseif abs(alpha(3)) <= 25
        ds3 = 1/BE * (1 - exp(alpha(3))/(1 + exp(alpha(3))));
    end
    %Chapman Parameters
    BC = 1 - 10 * (s(1)*ds1 + s(2)*ds2 + s(3)*ds3) / (s(1) + s(2) + s(3));
    z = (h - 100) / 10;
    %Electron Density
    Nbot = (s(1) + s(2) + s(3)) * exp(1 - BC*z - exp(-z)) * 1e11;
end

end